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ABSTRACT 

In a previous paper we formulated the problem of the formation and evolution of frag- 
ments (or cores) in magnetically-supported, self-gravitating molecular clouds in axisymmet- 
ric geometry, accounting for the effects of ambipolar diffusion and Ohmic dissipation, grain 
chemistry and dynamics, and radiative transfer Here we present results of star formation simu- 
lations that accurately track the evolution of a protostellar fragment over eleven orders of mag- 
nitude in density (from 300 cm^-^ to « 10^'* cm^'^), i.e., from the early ambipolar-diffusion- 
initiated fragmentation phase, through the magnetically supercritical, dynamical-contraction 
phase and the subsequent magnetic decoupling stage, to the formation of a protostellar core 
in near hydrostatic equilibrium. As found by Fi edler & MouschoviasI (1 1993b . gravitationally- 
driven ambipolar diffusion leads to the formation and subsequent dynamic contraction of a 
magnetically supercritical core. Moreover, we find that ambipolar diffusion, not Ohmic dis- 
sipation, is responsible for decoupling all the species except the electrons from the magnetic 
field, by a density « 3 x 10^^ cm"''. Magnetic decoupling precedes the formation of a central 
stellar object and ultimately gives rise to a conce ntration of magnetic flux (a ' magnetic wall') 
outside the hydrostatic core — as also found bv lTassis & MouschoviasI (l200 5a b) through a 
different approach. At approximately the same density at which Ohmic dissipation becomes 
more important than ambipolar diffusion (> 7 x 10^-^ cm^'^), the grains carry most of the 
electric charge as well as the electric current. The prestellar core remains disclike down to 
radii ~ 10 AU, inside which thermal pressure becomes important. The magnetic flux prob- 
lem of star formation is resolved for at least strongly magnetic newborn stars by this stage of 
the evolution, i.e., by a central density « 10^* cm^'^. The hydrostatic core has radius « 2 AU, 
density w 10^'' cm~'', temperature « 300 K, magnetic field strength « 0.2 G, magnetic flux 
« 5 X 10^* Wb, luminosity 10~^ Lq, and mass 10~2 M©. 

Key words: ISM: clouds — magnetic fields — MHD — stars: formation — radiative transfer 
— dust, extinction 



1 INTRODUCTION 



In a previous paper jKunz & MouschoviasI [2OO9I. hereafter Paper 
I) we formulated the problem of the formation and subsequent 
evolution of fragments (or cores) in magnetically-supported, self- 
gravitating molecular clouds in two spatial dimensions. The for- 
mulation allowed for the detailed study of both ambipolar diffusion 
and Ohmic dissipation, grain chemistry and dynamics, and radia- 
tion. Here we present results from the first star formation simula- 
tion to rigorously and simultaneously take into account both non- 
ideal magnetohydrodynamic (MHD) effects and radiative transfer. 

The thermodynamic evolution of protostellar cores in the 
absence of magneti c fields has been s tudied numeri c ally in 
spherical geometry ( iBodenheimeJ 1 19681 : iLarsorJ 1 19691 . Il972bl : 
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[ Appenzeller & Tschamutej ^ 1 19751 : lYorke & Krugell 1 19771: lYorkel 
1979: Win kler & Newm an 1980a b: Stabler. Shu & Taanilll980allbl: 
iMasunaga. Mivama & Inutsukal |1998 : Lesaffre et al.* ^IWy), in 
axisy n imetric geometry^ ( |Larson 1972al: iBlack & Bodenheiniei 



1976: Tscharnuter 1975, 1978: Boss 1984: Bodenheimer et al 



1990; Yorke,JBodenheimer & Laughlin 1993, 1995 ), and in three- 
dimensional geometrv JBossI [19 86. 1988. Jl993l : IBoss & Mvliiiij 
I1995I : lBate|[l998l : IWhitehouse & Batei.200a) . using widely vary- 
ing initial conditions, approximations, and numerical techniques. 
Despite these differences, a consensus has emerged regarding the 
qualitative evolution of hydrodynamic, nonisothermal prestellar 
objects that i s in broad agreement with the pioneering work of 
iLarsonl ( IT969I) . Such objects have been shown to remain isother- 
mal until central densities in the approximate range 2 x 10* cm~'^ 
- 2 X 10^" cm~'^ (depending on the simulation) are reached. The 
subsequent rise in temperature ultimately leads to the formation 
of a hydrostatic core at a temperature ^ 100 K and central den- 
sity ~ 2 X lO^'^ cm~^. The initial mass and radius of the hydro- 
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static core are typically found to be ~ 0.005 and ~ 4 AU, 
respectively. A thermal shock forms at the core boundary as mass 
is accreted onto the hydrostatic core, while small-amplitude oscil- 
lations occur about equilibrium. Finally, a temperature of 1000 K 
is reached at a central density ~ 2 x 10^^ cm~'^. Soon thereafter, 
molecular hydrogen dissociates (at 2000 K), the ratio of specific 
heats 7 dips below the critical value 4/3, and material at the centre 
of the core becomes unstable and collapses dynamically. Despite 
the ensuing highly dynamic evolution, the temperature rises only 
slowly since most of the gravitational energy goes into molecular 
dissociation. Once molecular hydrogen is almost all dissociated, 
7 rises above 4/3 and the thermal pressure increases rapidly, de- 
celerating and ultimately halting the collapse at the centre. A sec- 
ond (stellar) core forms, accompanied by another small rebound 
and subsequent radial pulsations. The mass and size of this core 
are « 0.0015 and « 1.3 R©, respectively. Its central density 
« 0.02 g cm"^ and temperature f» 20, 000 K. 

The addition of magnetic fields to the picture has brought a 
myriad of interesting effects (e.g. , see reviews by iMouschoviaj 
1 1981 1 1994 ICiolek &Basi3l2000l) . Aside from the support that 
magnetic fields provide against self-gravity, the inclusion of mag- 
netic fields has led to a qualitatively different picture of molecu- 
lar cloud fragmentation and evolution than earlier hydrodynamic 
calculations had suggested. In an initially magnetically subcritical 
cloud, gravitational infall occurs only as rapidly as allowed by am- 
bipolar diffusion, the relative drift of neutral and charged particles. 
Ambipolar diffusion allows thermally supercritical fragments to 
separate out, in which the neutral particles fall in toward local grav- 
ity centres, impeded by collisions with the charged particles, which 
are essentially 'held in place' by magnetic forces. Magnetic braking 
is very effective during this phase of contraction and is responsible 
for reducing the angular momenta of cloud cores to their observed 
low values and for resolving the angular momentum problem dur- 
ing the early, isothermal stage of contraction. The central mass-to- 
flux ratio of each forming fragment eventually increases to the crit- 
ical value for collapse (typically when the central density reaches 
values in the range ~ 10* cm~"^ - 10^ cm"''). The now thermally 
and magnetically supercritical fragments begin to collapse more 
rapidly than their surroundings and they are simply referred to as 
supercritical fragments or cores. Each core evolves dynamically 
(though slower than free fall) under near flux freezing, until the 
resurrection of ambipolar diffusion causes magnetic decoupling] 
to set in at a density ~ 10^" cm~"^. Magnetic decoupling occurs 
over a few orders of magnitude in central density enhancement and 
precedes the formation of a central stellar object. The isothermal 
evolution of a magnetic molecular cloud has been followed in detail 
up to c entral densities of 2 x 10^^ cm~'^ bv lOesch & MouschoviasI 
( I2OOII) , who find that the magnetic field strength asymptotically ap- 
proaches ~ 0.1 G in the innermost ~ 20 AU of the c l oud. 

The nu merical simulations of lBossI ( Il99l 1 19991 , l2002l 12001 
I2OOII2OO9I) followed the thermodynamic evolution of a magnetic 
protostellar core with a reasonable degree of accuracy, via the Ed- 
dington approximation, while treating magnetic field effects and 
ambipolar diffusion crudely through various approximations and 
parametrisations based on prev ious isothermal MHD calculations. 
iTassis & MouschoviasI j2007al lbll3) took the opposite approach by 



^ Complete magnetic decoupling refers to conditions such that the mag- 
netic field has no significant effect on the dynamics of the neutral matter 
and the motion o f the neutral matter no longer affects the magnetic field — 
see footnote 3 of lDesch & MouschoviasI i200ll) . 



incorporating a detailed treatment of nonideal MHD and chemistry 
into their numerical simulations, while approximating the ther- 
modynamic evolution via a piec ewise adiabatic equation of state. 
More recently, Com mercon et aL ( 2009) followed the evolution of a 

1 M0 dense core in the presence of magnetic fields, including flux- 
limited diffusion radiative transfer but neglecting both nonideal- 
MHD effects and chemistry. 

The present paper follows the detailed nonideal MHD and 
grain physics and chemistry of Tassis & M ouschovias ( 2007a, b, 3), 
but also adds radiative transfer, which is essential for making accu- 
rate predictions that can be compared with observations pertaining 
to the opaque phase of star formation. In Section|2]we summarise 
the formulation of the problem and method of solution (detailed 
in Paper I), and we provide the initial and boundary conditions of 
the numerical simulation. Section [3] presents the results. Contact 
is made with both observations and prior theoretical work where 
appropriate. Section|4|summarises the results and predictions, and 
discusses their limitations. 

2 PHYSICAL MODEL 

A detailed discussion of the physical processes included in the cal- 
culations was presented in Paper I. In summary, we consider a two- 
dimensional, nonrotating, nonisothermal model molecular cloud, 
whose axis of symmetry is aligned with the z-axis of a cylindrical- 
polar coordinate system (r, (j>, z). The cloud is initially threaded by 
a uniform magnetic field oriented along the symmetry axis. The 
abundances of all species (except the neutrals) are determined from 
an extensive chemical-equilibrium network that includes UV radia- 
tion, cosmic rays, radioactive decays, thermal ionisation, dissocia- 
tive and radiative recombination, atomic and molecular ion charge 
transfer, electron and ion attachment onto grains, and charge trans- 
fer by grain-grain collisions. While the code described in Paper I 
can handle an MRN grain size distribution, we present here only the 
simpler case of a uniform grain size distribution. Results of numeri- 
cal calculations that employ an MRN distribution will be presented 
as part of a parameter study in a future publication. 

The model cloud is evolved by using the six-fluid RMHD 
equations (Equations 72a-g in Paper I) from a nonequilibrium, uni- 
form 'reference' state, characterised by a number density of neu- 
trals nn,ref = 300 cm~"^ (corresponding to a neutral mass density 
Pn.rof ^ 1.17 X 10"^^ gcm"^), temperature Trof = 10 K (cor- 
responding to an isothermal sound speed Crcf = 0.188 km s"^), 
and magnetic field strength B^cf ~ 15 fiG (corresponding to an 
Alfven speed va ~ 1.24 km s~^). The part of the cloud whose 
evolution is followed numerically has radius R and half-thickness 
Z, both equal to 0.75 pc, implying a total mass of 45.5 Mq. The 
central mass-to-flux ratio of the reference state in units of the criti- 
cal central value is then 

(rfM/rf$B)c ^ 2Z/0n,rcf/grof 

~ (dAf/dcE.B)c,cr " (3/2)(63G)-i/2 

where (dM/d$B)c,cr = (3/2)(63G)~^' ^^ is the central critical 
value f or collapse originally determined bv lMouschovias & Spitzej 
(fl97^ . 

If the magnetic field were frozen in the matter, the 'reference' 
state would relax along field lines and oscillate about an equilib- 
rium state (denoted by a subscript '0') in which gravity is balanced 
by thermal-pressure forces along field lines and mainly by magnetic 
forces perpendicular to field lines. An important dimensionless pa- 
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rameter in this equilibrium state is tlie ratio of tiie magnetic and 
thermal pressures. 
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In equilibrium (see iMouschoviasll 1 99 Ibl) , these two dimensionless 
free parameters satisfy the relation: 



Qo/io ~ 0.71. 



(3) 



The physical meaning of this equation is that, since the same force 
(gravity) is balanced by the thermal-pressure force along field lines 
and by the magnetic force perpendicular to the field lines, the 
thermal-pressure and magnetic forces must also be comparable in 
magnitude. A point related to this equation cannot be overempha- 
sised: for subcritical central flux tubes (/io < 1), the local mag- 
netic pressure exceeds the thermal and gravitational pressures. This 
does not imply, however, that such flux tubes will expand. It sim- 
ply means that the magnetic field in the core is comparable to that 
in the envelope, and is confined by the massive cloud envelope, 
not the gravity of the core. Virial-theorem based arguments miss 
this essential point, leading to the incorrect conclusion that mag- 
netically subcritical fragments (or clouds) cannot exist. After all, 
iMouschov ias ( 1976a b) calculated rigorously the equili br ia of s uch 
objects, as did lTomisaka. Ikeuchi & Nakamural ( ll988allbllT98^ . 

Since /x^cf = l-^o (i.e., contraction along field lines does not 
alter a fragment's mass-to-flux ratio), the central density in the 
equilibrium state will be a factor ~ lAaicifJ-r^t ~ 7 greater than 
that of the reference state. As in Fiedler & Mouschovias ( 199:1), 
we have found little quantitative difference (for typical parame- 
ters) between a run in which ambipolar diffusion operates from 
the outset and one in which the cloud is allowed to reach equilib- 
rium under flux-freezing before ambipolar diffusion is switched on. 
The physical reason for this is that the relaxation time along field 
lines is short compared to the ambipolar diffusion timescale, par- 
ticularly for the low densities characteristic of the relaxation phase 



{Un < 10'' cm" 



see below), during which ionisation due to 



UV radiation is important. 

Once the heat generated by released gravitational energy dur- 
ing core contraction is unable to escape freely (at a central number 
density riopq — lO'^ cm~^), we use radiative transfer to determine 
the thermal evolution of the core and its effect on the dynamics. For 
this, we employ the grey flux-limited diffusion (FLD) approxima- 
tion. Realistic temperature-dependent grain opacities are used that 
account for a variety of grain compositions (see Paper I, fig. 2). 

Due to the assumed symmetry of the problem, we need only 
follow the evolution of the region bounded by ^ r ^ i? and 
^ z !^ Z. Boundary conditions on the density and velocity are 
chosen so that no mass, magnetic flux, or thermal energy crosses 
any boundar y, as employed and tested bv iFiedler & Mouschovias! 
( Il99lll993h . The magnetic field lines are constrained to cross the 
midplane (z = 0) and upper boundary (z — Z) normally, and to 
have no radial component at the axis (r = 0) and outer boundary 
(r = R). In order to prevent unrealistic densities near the upper 
boundary of the computational domain (z = Z), a floor is imposed 
on the neutral number density at rin, floor ~ 1 cm~'^. In addition, a 
floor of 10 K is imposed on the gas and radiation temperatures. 

The number of computational cells is fixed (80 x 80), and 
their positions are spaced logarithmically, so that the spacing be- 
tween the ith and ith + 1 cells is a number greater than the spacing 
between the ith — 1 and ith cells. This ratio is kept spatially uni- 



form (with the same value in both the r and z directions) but is 
allowed to vary in time. Its value is such that the innermost cell 
is constrained to always have a width AT,cr/5 - A T.cr/7, where 
Ax.cr = IACts is the critical thermal lengthscale jMouschoviasI 
Il991al) and Tfl = (37r/32G/9n.c)^''^ is the central spherical free- 
fall time. (The quantity pn,c is the neutral mass density at the cloud 
centre, and the local isothermal sound speed C in this calculation is 
a function of time.) The critical thermal lengthscale is the smallest 
scale on which there can be spatial structure in the density without 
thermal-pressure forces smoothing it out. At the end of the simula- 
tion, the minimum cell size is ~ 0.2 AU x 0.2 AU. 



3 RESULTS 

3.1 Overall Evolution 

In Fig. [T] we show (a) the time evolution of the central number 
density of neutrals, nn,c, and (b) the central mass-to-flux ratio, 
(dm/d<l>B)c, (normalised to the central critical value for collapse) 
as a function of Un.c- (Note that our normalisation is greater than 
that of Tassis & MouschoviasI ( 12007 al lbll3). which referred to a uni- 
form, thin disc, by a factor « 2.) There are three distinct phases of 
evolution: relaxation, quasistatic, and dynamic. First, the cloud re- 
laxes along magnetic field lines from its uniform reference state to a 
quasi-equilibrium state whose central density rin.c — 2280 cm"^. 
Ambipolar diffusion is negligible during this phase, since it is char- 
acterised by relatively low densities (< 10"^ cm~^) and therefore 
relatively large degrees of ionisation (~ 10"'' - lO"'^) due to UV 
radiation. Hence, the mass-to-flux ratio remains roughly constant. 
Then ambipolar diffusion sets in and the fragment contracts qua- 
sistatically (i.e., with negligible acceleration) under its own self- 
gravity until a supercritical core forms at nn,c — 1.1 x 10^ cm~'^ 
and time « 9 Myr (marked in the figure by a 'star'). The subse- 
quent evolution is dynamic, although significantly slower than free 
fall. The central mass-to-flux ratio asymptotes to roughly twice its 
critical value until ambipolar diffusion is reawakened at a central 
density nn,c ~ 10^^ cm~^. The mass-to-flux ratio then increases 
dramatically, reaching ~ 80 times the central critical value for col- 
lapse by the end of the calculation (at a density ~ 10^* cm""'). 
During the entire supercritical phase of evolution, the central den- 
sity increases by ten orders of magnitude in 2.38 Myr. 

Note that soon after the core becomes magnetically super- 
critical (at a central density 1.1 x lO'* cm"''), it will also be- 
come observable as a well-defined ammonia core. Therefore, Zee- 
man observations of ammonia cores are predicted to preferentially 
measure slightly supercritical mass-to-flux ratios. This is in excel- 
lent agreement with observations (e.g., by Crutcher 1999). Also, 
the detectable (using ammonia as a tracer) lifetime of the core 
formed in this simulation is ~ 1 Myr, which is also in excellent 
ag reement with the ^ 1 My r lifetimes of starless cores adv ocated 
by 'Lee & Myers ('1999^ and ljiiina. Myers & Adain3 ( Il999l) . (See 
fPassis & Mouschovias 2004] for further discussion of this point and 
its importance in interpreting starless core statistics.) 

The overall spatial and temporal evolution of the cloud is 
shown in Figures |2] and [3] In Fig. |2] the density (thin solid lines), 
magnetic field (thick solid lines), and velocity field (arrows) are 
displayed in each of the eight frames, which show portions of 
the cloud at different times, when the central density nn,c ~ 
3 x (10^ lO'', 10^ 10^ 10^ 10^\ lO", lO^^) cm-^* (leftto 
right, top to bottom). Fig.[3]shows the innermost 0.01% of the cloud 
at the end of the simulation (when nn,c — 10^^ cm"''). In Figures|2] 
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Figure 1. (a) Central number density of neutrals, 7in,c, as a function of 
time, (b) Central mass-to-flux ratio (normalised to the central critical value 
for collapse) as a function of rin.c- The 'star' marks the time at which a 
supercritical core forms. 



and[3^, the velocity vectors are normalised to the maximum veloc- 
ity in each frame (given below) and every third isodensity contour 
denotes a change in density by a factor of 10. In Fig.[3}5, the den- 
sity (dashed lines) and temperature (solid lines) are displayed, with 
every tenth isotherm denoting a change in temperature by a factor 
of 10. 

In the first row of Fig. [2] we show the entire computational 
domain (Q ^ r ^ R, Q ^ z ^ Z) ?i the three instances at which 
the central density is nn,c = 3 x (10^, 10'', 10^) cm"^. In the 
first frame, the cloud has settled into a quasi-equilibrium state 
and ambipolar diffusion has already set in. The magnetic field 
lines remain relatively undeformed, while the neutrals contract 
through them under their self-gravity - their maximum radial 
velocity at this moment is 0.012 kms^'). For r > 0.5 pc and 
z > 0.2 pc, there is a small outward motion as the cloud rebounds 
from overshooting its equilibrium state. In the second frame, 
the central magnetic flux tubes of the cloud have just become 
magnetically supercritical and begun to bend into an hourglass 



morphology. Such a field morphology has be en observed via 
dust p o larimetry by, e.g., Schleuning ( 1998), Hildebrand et al.l 
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19990, iGirart. Crutch er & Rao ( 1999), Schleuning et af ( 200( 



iLaietalJ ilOOj). 



Matthews & Wilson (2002), Houde et al 

(l2004l).ICortes & Crutcheil l2006l), lOirart Rao & Marr one (200^ 



IVaillancourt et al.l ( l2008l) 7 lTang et al.l ( l2009h . and lKirbv, 1 .200% , 



The maximum radial velocity has increased to 0.033 km s" 
The outermost envelope of the cloud remains magnetically 
supported, whence the radial velocities there are very small. A 
shock has formed at z ~ 0.3 pc (maximum vertical velocity 
0.39 kms~' ~ 2C) because of the rapid contraction along field 
lines. In the third frame, the magnetically supercritical core is well 
into its dynamic (though slower than free-fall) stage of evolution 
and the maximum radial velocity has reached ~ 0.1 km s~'. 

By this point, the qualitative features of the outer envelope 
of the cloud, which is magnetically supported, have been decided 
and hardly change at subsequent times. The result is a highly non- 
homologous collapse. We therefore focus on progressively smaller 
radii and track the local evolution of the magnetically supercritical 
core. In the first frame of the second row, we show the innermost 
10% of the computational region when Un.c = 3 x lO*" cm~^. The 
maximum radial (vertical) velocity at this moment is 0.18 km 
(0.47 km s~'). It is clear that the core remains disclike down to 
the smallest scales, even in the absence of rotation, a consequence 
of the presence of magnetic forces. In the second frame of this 
row, we show the innermost 1% of the computational region when 
"-n,c = 3 X lO'"* cm^"^. The maximum radial (vertical) velocity 
at this moment is 0.38 kms^^ (0.80 kms^^). The similarities 
between these two figures suggest that the supercritical phase of 
core con traction is ne a rly self-similar, a fact that has been exp loited 
by, e.g., iBa su' ('1997), 'Contopou los, Ciolek & Konigll ( Il998l) , and 
iKrasnoDolskv & Konigl (2002) , who derived semi-analytic self- 
similar solutions for this stage of contraction that reproduce the 
main qualitative features found in detailed simulations. 

The evolution departs from being nearly self-similar at higher 
densities, once the gas begins to decouple from the magnetic field 
and a transition from disclike to spherical geometry takes place. 
For the third row of frames (rin.c = 3 x lO" , 10^^ 10^^ cm"^), 
we zoom in on the innermost 0.01% of the computational region 
(~ 15.5 AU X 15.5 AU) in order to highlight the formation of the 
hydrostatic core. In the first frame (nn,c = 3 x 10^^ cm^^), the 
maximum radial (vertical) velocity is 0.52 km (1.03 km s~'). 
This is the last frame shown in which the innermost r < 4 AU 
of the cloud exhibits a disclike geometry. In the second frame 
(nn,c = 3 X 10^^ cm~^), this region has begun to assume a spher- 
ical geometry due to the increasing relative importance of thermal 
pressure. The magnetic field has essentially decoupled from the 
matter and the field lines are straight inside the ~ 3 x 10^' cm""^ 
isodensity contour. The radial and vertical velocities attain their 
maxima (0.64 km and 1.14 km s~^, respectively) just outside 
the boundary of the newly-formed hydrostatic core, where shocks 
occur in both radial and vertical directions. These shocks are a re- 
sult of the rapid deceleration of matter as it slams into the hydro- 
static core boundary. In the third frame (nn,c = 3 x 10^^ cm""^), 
the hydrostatic core is clearly visible and exhibits a spherical ge- 
ometry inside ~ 2 AU. The magnetic field lines continue to re- 
main straight within the ~ 3 x 10^' cm""^ isodensity contour. The 
maximum radial (vertical) velocity in the frame is 1.13 kms^' 
(1.27 kms~'). A closer inspection reveals a small outward mo- 
tion with speed of ~ 0.1 km s"'^ along the vertical symmetry axis, 
a result of infalling gas rebounding after overshooting hydrostatic 
quasi-equilibrium. Throughout the evolution, the magnitude of the 
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Figure 2. Overall spatial and temporal evolution of the model molecular cloud. The density (thin solid lines), magnetic field (thick solid lines), and 
velocity field (arrows) are given in each of the eight frames, which show portions of the cloud at the times when the central density nn,c = 3 x 
(10^, 10*, 10^, 10®, 10^, 10^1, 10^^, 10^^) (left to right, top to bottom). The velocity vectors are normalised to the maximum velocity in 

each frame (see text). Every third isodensity contour denotes a change in density by a factor of 10. 
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Figure 3. Overall structure of the innermost 0.01% of the model molecular 
cloud at the end of simulation (nn.c = 10^* cm~'^). (a) Density (thin 
solid lines), magnetic field (thick solid lines), and velocity field (arrows). 
The velocity vectors are normalised to the maximum velocity in the frame 
(1.44kms-l); every third isodensity contour indicates a change in density 
by a factor of 10. (b) Density (dashed lines) and temperature (solid lines). 
Every tenth isotherm denotes a change in temperature by a factor of 10. The 
maximum temperature is 266 K. 



radial component of the magnetic field never becomes greater than 
the magnitude of the vertical component. Magnetic pinching forces 
remain small and magnetic reconnection does not occur. 

Fig. [3^ shows the density, magnetic field, and velocity field 
in the innermost 0.01% of the computational region at the end 
of the simulation (when Wn.c — 10^* cm""^). The maximum ra- 
dial (vertical) velocity in the frame is 1.44 km (1.24 km s~^). 
In Fig. [3j), we highlight the temperature structure by overplotting 
isodensity contours (dashed lines) and isotherms (solid lines). The 
central temperature at this point is 266 K. Isotherms and isodensity 
contours coincide only within a region ^= (r^ + z^)^''^ < 2 AU. 
Outside this region, the temperature structure is appreciably more 
spherical than the density structure. 



3.2 Evolution of Central Quantities 

Fig. exhibits the central gas (solid line) and radiation (dashed 
line) temperatures, Tc, as a function of the central number den- 
sity of neutrals. The gas remains isothermal at 10 K until a density 
riopq — 10^ cm^"^, after which the heat released due gravitational 
contraction is unable to escape freely. The central gas tempera- 
ture then rises with density at a rate 91nTc/91nnn,c = 7eff — 1, 
where 7cft is the effective adiabatic index denoted by the solid line 
in Fig. I4J3. The dotted line in Fig. |4j3 shows this index if the gas 
were to evolve adiabatically beyond a central density 10^^ cm^^. 
The gas never evolves adiabatically during the simulation. There is 
some indication that adiabaticity will set in at the cloud centre by a 
central density ~ 10^^ cm"'', after which 7cft — 7/5. 

A central temperature of 100 K is reached when the central 
density is ~ 6 x 10^^ cm "''. This is a si gnificantly higher den- 
sity than that found by, e.g.. lLarsoril ( Il969t) . There are two princi- 
pal reasons for this difference. First, a disclike, rather than spher- 
ical, geometry allows the radiation to escape more easily, since 
the radiation is not isotropically confined by a spherical opacity 
distribution. Second, the early redistribution of mass by ambipo- 
lar diffusion is responsible for lowering the dust-to-gas ratio in 
the central flux tubes of the cloud from its usual interstellar value 
jCiolek & Mouschovias 1994. ; also discussion below). The fewer 
dust grains there result in a lowered opacity, and consequently al- 
low radiation to escape more easily. This effect is complicated 
somewhat, however, as our op acities are t emperature-dependent 
(see fig. 2 of Paper I), whereas LarsonI ( Il96 9) assumed a constant 
value xr = 15 cm^ of dust (with a constant dust-to-gas ratio 
of 0.01). 

Fig. |4j; exhibits the central magnetic field strength, Be, as 
a function of central density. During the ambipolar-diffusion- 
controlled evolution (until the formation of a supercritical core, 
marked by the 'star'), Be increases by only ~ 47%, from 15 ^lG 
to 22 fj,G. A phase of near flux-freezing follows, during which Be 
resembles a power-law Be oc n(^.e. with k ~ 0.47 (see Fig.|4ll), 
as found bv iFiedler & MouschoviasI ( Il993h . This is in excellent 
agreement with the value k = 0.47 ±0.08 inferred from Zee- 
man observations of protostellar cores ( ICrutchei|[T99Sl) . For a thin 
disc in the flux-freezing limit. Be oc nU^cTc^'^ , so that k = 7off/2 
( lMouschoviaslll976allbril991bh . Therefore, any nonisothermal evo- 
lution (i.e., 7eff > 1) will result in k > 1/2 (in the flux-freezing 
limit). This is evident briefly in Fig.|4ll for the central density range 
"-n,e ~ 10^ - lO'^^ cm""'. This trend is halted once the magnetic 
field begins to decouple from the matter and k decreases. 

Note that the rapid, large oscillations in k seen in fig. If of 
iTassis & MouschoviasI ( l2007bh after magnetic decoupling sets in 
are absent in our Fig.|4|l. Those oscillations were due to the abrupt 
and (relatively) rapid adiabatic temperature increase in the core, 
which forced the contracting core to overshoot and subsequently 
oscillate about hydrostatic equilibrium. In the situation presented 
here, the much more gradual increase in the central temperature 
affords the core enough time to maintain hydrostatic force balance 
as it contracts. A close inspection of Fig. |4jl does indeed reveal 
small physical oscillations about a hydrostatic equilibrium. 

The evolution of the central (gas) temperature assists in under- 
standing the evolution of the quantity 91n En.e/91n nn,c (shown 
in Fig. Is}, which is a quantitative indicator of the core geome- 
try. For a disc, in which thermal-pressure forces balance gravity 
along magnetic field lines, or a sphere, in which thermal-pressure 
forces balance (spherical) gravity radially, 91nEn,c/91nnn,c = 
7off/2 (whose value is represented by the dotted line in the figure), 
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Figure 4. Evolution of central quantities as functions of the central number density of neutrals, n^^c' (a) Gas (solid line) and radiation (dashed line) tempera- 
tures, (b) Exponent 7(,ff in the relation Tc oc rtn°J^ ^ (solid line). The dotted line shows the exponent if the gas were to evolve adiabatically beyond a central 
density 10^^ cm~'^. (c) Magnetic field strength, (d) Exponent k in the relation Be oc nJJ ^ (solid line). The dotted line shows the exponent if the magnetic 
field were to evolve under strict flux-freezing and the gas were to behave adiabatically beyond a central density 10^^ cm^'^. The 'star' marks the time at 
which a supercritical core forms. 



whereas for a sphere of constant mass, d In En,c/9 In nn,c = 2/3. 
After the rapid contraction along field lines, 91nEn,c/91nnn,c 
approaches and oscillates about 1/2; the cloud has settled into 
a disclike quasi-equilibrium. Even during the magnetically su- 
percritical phase of dynamical contraction, the geometry remains 
disclike (as indicated by the near coincidence of the solid and dot- 
ted curves in the figure). This trend holds until central densities 
3 X 10^^ cm~'^ are attained, after which departures from disc 
geometry become pronounced and, by the end of the run, the for- 
mation of an almost spherical core is evident. 

Fig. [6] shows the evolution of quantities related to the chem- 
istry as functions of central density: (a) the ionisation rate at the 
centre of the cloud, ("c, due to all processes (UV radiation, cos- 
mic rays, and radioactive decays); and (b) the central abundances 
of different species, Xb,c = Ws.c/nn.c (where s = e, i, g_, go, 
g+). For Un.c 1800 cm~^, UV radiation dominates the ion- 
isation rate (see Paper I, § 4.1). Beyond this density, the ionisa- 



tion in the core is mainly due to cosmic rays. Its value is constant, 
( = 5 X 10"" s"\ for densities up to 3 x lO" cm"^, be- 
yond which the column density exceeds ~ 100 g cm^^ and cos- 
mic rays are appreciably attenuated. The ionisation rate then de- 
creases monotonically, until it is dominated by radioactive decays 
(mainly of ''"K) above a central density ~ 10^'^ cm~"^ and reaches 
a constant value = 6.9 x 10"^'^ s~^. Thermal ionisation does 
not become important for the density range investigated here, since 
the temperature never reaches ~ 10'^ K. A simple extrapolation of 
temperature based on our results indicates that such a temperature 
will not be reached until a central density ~ 3 x 10^^ cm""^. It is 
therefore suggested that magnetic recoupling will not occur until 
central densities of at least several x 10^^ cm^^ are attained. 

The central species abundances are strongly affected not only 
by changes in the ionisation rate, but also by microscopic interac- 
tions between the different species and by the macroscopic dynam- 
ics of the cloud. [For central neutral densities < 1800 cm~"^, the 
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10" IC" 10^ 10'° 10^^ 10^" 
^n,c (cn-r^) 

Figure 5. Evolution of the exponent In Snjc/c* lnnn,c (solid line) as a 
function of the central number density of neutrals, nn,c- The dotted line 
shows the exponent if the geometry were to remain disclike throughout the 
evolution. The 'star' marks the time at which a supercritical core forms. 



term added to the ion density (as in lFiedler & Mouschoviasiri992l 
Il993l) . which mimics the effect of UV ionisation at low densities 
in cloud envelopes, is responsible for the behaviour of the ion and 
electron relative abundances and is not meant to be a model or pre- 
diction for the exact behaviour of ni and in cloud envelopes. It 
is used mostly for computational convenience.] 

At densities greater than ~ 2 x 10^ cm^^ ionisation is pri- 
marily due to cosmic r ays, a nd Xi,c behave as found earlier by 
ICiolek & Mouschovias! ( Il994l) . It cannot be overemphasised that 
the 'canonical' relation Xi,e oc n^"^ with k = const — 1/2 
is never established. In fact, during the subcritical phase of the 
evolution (i.e., nn,c ^ 10* cm~^), k is always greater than 1/2, 
an effect that is entirely due to ambipo lar diffusion. As explained 
bv lCiolek & Mouschovial ( Il994 [ 1998h . during this phase of evo- 
lution charged and neutral dust grains are well atta ched to mag- 
netic f ield li nes, which are essentially 'held in place' ( iMouschoviasI 
Il978lll979l) , and are 'left behind' by the inwardly diffusing neu- 
trals. This results in a decrease in the central dust-to-gas ratio by 
a factor roughly equal to the initial (dimensionless) central mass- 
to-flux ratio. One consequence is a significantly reduced inelastic 
capture of ions and electrons onto grains and therefore a greater 
number of gas-phase ions and electrons than predicted by calcula- 
tions that assume a constant dust-to-gas ratio. 

Once a supercritical core forms, the grain relative abundances 
are frozen at their values until a central density of ~ 10^ cm""' 
is reached. Gas-phase ions and electrons are then quickly adsorbed 
onto grain surfaces and the grains become the main charge car- 
riers (when rin.c ~ 10^" cm^"^). However, most grains remain 
neutral due to inadequate numbers of ions and electrons, whose 
abundances are determined by balancing their rates of production 
and their collision rates with neutral grains. Because these collision 
rates are inversely proportional to the square root of the ion/electron 
mass, the ions are more abundant than the lighter electrons by a 
factor ~ (rrii/mo)"''^ — 210 (different sticking probabilities for 
ions and electrons reduce this value to ~ 125). Since x^^ ap- 
proaches a constant, the densities of ions and electrons approach 
constant values, so their abundances relative to neutrals asymptote 




10" 10** lO" 10^° 10'^ 10^" 
■^n,c (cm ) 

Figure 6. Evolution of central quantities as functions of the central number 
density of neutrals, nn,c: (a) Ionisation rate due to all processes (UV radia- 
tion, cosmic rays, and radioactive decays), (b) Abundances of species. The 
'star' marks the time at which a supercritical core forms. 

towards a " dependence. Once cosmic rays become appreciably 
attenuated, however, this asymptotic behaviour is interrupted and 
the abundances of electrons and ions decrease much more rapidly. 
An ionisation floor is established by the radionuclide at den- 
sities ~ 10"^'^ cm^"^, after which the electron and ion abundances 
once again scale roughly as n^" . The abundances of charged grains 

— 1/2 

are proportional to Wn , for the same reasons that applied to the 
abundances of electrons and ions at low densities. 

Central quantities related to magnetic diffusion are shown in 
Fig. |7] as functions of central density. The magnetic diffusion co- 
efficients related to ambipolar diffusion (AD) and Ohmic dissipa- 
tion (OD) are shown in Fig. [7^. As we are concerned here with 
the evolution of only the poloidal component of the magnetic field 
(recall that rotation has been justifiably ignored), the Hall effect 
plays no role in the dynamics of the cloud. Magnetic-flux re- 
distribution is due to ambipolar diffusion until a central density 
~ 7 X lO"^ cm~"^ is reached, after which Ohmic dissipation is the 
primary magnetic diffusion mechanism. Therefore, ambipolar dif- 



Nonisothermal magnetic star formation 9 




Figure 7. Evolution of central quantities as functions of the central number density of neutrals, nn,c- (a) Resistivities due to ambipolar diffusion (AD) and 
Ohmic dissipation (OD). (b) Parallel (||), perpendicular and Hall (H) magnetic conductivities, (c) Fraction of the parallel magnetic conductivity carried 
by each species, (d) Fraction of the perpendicular magnetic conductivity carried by each species. The 'star' marks the time at which a supercritical core forms. 



fusion still dominates Ohmic dissipation as a flux reduction mecha- 
nism when the grains become the primary charge carri ers, contrary 
to the expectations of iNakano & Umebavashil ( Il986al lbl). Equality 
of the ambipolar-diffusion and Ohmic-dissipation diffusion rates 
occurs once the central mass-to-flux ratio becomes equal to ~ 16 
times its critical value and the grains become the dominant current 
carriers (see below). 

In Fig.|7j5, we show the parallel (||), perpendicular (±), and 
Hall (H) conductivities. Figures|7}; and|7}l show the fraction of the 
parallel and perpendicular conductivities, respectively, carried by 
each species. (Absolute values have been taken; the cusps in the 
positive grain curves are due to the conductivity passing through 
zero as it changes sign.) The parallel conductivity is primarily due 
to electrons until Ohmic dissipation becomes the primary magnetic 
diffusion mechanism, after which the contribution of the charged 
grains to the parallel conductivity dominates that of the electrons. 
The perpendicular conductivity is mainly due to ions up to a central 
density ~ 3 x lO'^'^ cm~^, after which the charged grains dominate 
the ions' contribution. In other words, the grains begin to carry most 



of the current at roughly the same time that Ohmic dissipation be- 
comes important (see Section |3.3.5l l. Star formation and protostel- 
lar disc calculations that study the phase when Ohmic dissipation 
becomes an important magnetic diffusion mechanism must include 
not only the chemistry but also the magnetohydrodynamics of dust 
grains. 



3.3 Evolution of Spatial Structure of tlie Model Cloud 

Figures [8] - [18] show radial midplane (solid lines) and vertical 
symmetry-axis (dashed lines) profiles of physical quantities at thir- 
teen different times. Each curve corresponds to a time tj at which 
the central density has increased by a factor of lO-* relative to 
the initial uniform 'reference' state (i.e., the central density is 
~ 3x10^^-' cm~'' at time tj). The last curve (corresponding to ti2) 
is an exception; it refers to the time at which nn,c — 10^'' cm~^. 
These times are to = -4.00, ti = 3.763974, t2 = 10.022421, 
ts = 10.890328, t4 = 11.021924, = 11.049141, te = 
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Figure 8. Radial midplane (solid line) and vertical symmetry-axis (dashed line) profiles of (a) the number density of neutrals and (c) the gas temperature 
at thirteen different times (to = -4.00, ti = 3.763974, t2 = 10.022421, = 10.890328, t4 = 11.021924, is = 11.049141, te = 11.056294, 
t7 = 11.058372, ts = 11.059021, tg = 11.059230, tio = 11.059310, hi = 11.059346, ti2 = 11.059357 Myr). Also shown are the radial derivatives 
along the midplane of (b) the density profile, 91nnn/c'lnr, and (d) the (gas) temperature profile, dlnT/dlnr. The inner (outer) 'star' on a curve, present 
only after a supercritical core forms, marks the initial (final) radius of the supercritical core. 



11.056294, t7 = 11.058372, ts = 11.059021, tg = 11.059230, 
tio = 11.059310, til = 11.059346, and ti2 = 11.059357 Myr. 

3.3.1 Density and Temperature 

Fig.[8ji displays the number density of neutrals as a function of ra- 
dius (r) and height (z) at the 13 different times referred to above. 
We also show the logarithmic slope of the (midplane) density pro- 
file versus radius, 9 In Tin /9 In r, in Fig. [8]). A uniform central re- 
gion is maintained (in the radial direction) for r < AT,cr(t) by 
thermal-pressure forces. Outside this region, the density profile re- 
sembles a broken power-law. Inside the supercritical core, the slope 
has an average value ~ —2.2 for 10"* pc < r < lO"'^ pc. For 
smaller radii in the range 10"^ - 10"* pc 2 - 20 AU), the 
slope increases briefly due to a local increase in magnetic support 
(see Section [3.3. 3t . A transition to spherical geometry occurs at 
r ~ 4 AU when the central density is ~ 3 x 10*^ cm"''. The 



radial and vertical density profiles are nearly identical inside this 
radius, with a slope approaching ~ —10/3, corresponding to hy- 
drostatic force balance (with jcs = 7/5). Outside the magnetically 
supercritical core (~ 0.4 pc), there is an abrupt break in the slope 
because of significant magnetic support in the envelope. Very near 
the radial boundary (r = R = 0.75 pc) of the computational re- 
gion, the slope decreases again to satisfy the boundary condition 
that the matter has no radial velocity at this surface. 

Fig. [8}; displays similar profiles for the (gas) temperature. 
The logarithmic slope of the (midplane) temperature versus radius, 
91nT/91nr, is shown in Fig. [8}l. The molecular cloud and su- 
percritical core evolve isothermally for densities rin < 10^ cm""' 
and (cylindrical polar) radii r > 10"^ pc. Just inside the hy- 
drostatic core (r < 2 AU), the slope is closely approximated by 
(7off — l)(91nnn/91nr), whereas the temperature becomes con- 
stant very near the centre of the core. For 2 AU < r < 2000 AU, 
no single power-law can approximate the temperature profile. 
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Figure 9. Radial midplane (solid line) and vertical symmetry-axis (dashed line) profiles of (a) the 2-component of the magnetic field and (d) the ratio of the 
Alfven speed and the local isothermal sound speed at different times, as in Fig.[8] Also shown are (b) the radial derivative of the z-component of the magnetic 
field, d In /9 In r, along the midplane, and (c) radial profiles of the maximum (in z) strength of the r-component of the magnetic field. The inner (outer) 
'star' on a curve, present only after a supercritical core forms, marks the initial (final) radius of the supercritical core. 



3.3.2 Magnetic Field and Alfven Speed 

The z-component of the magnetic field is shown in Fig.|9ji. In com- 
mon with the density profile, there is an inner region [r < AT,cr(t)] 
in which the magnetic field is almost uniform. Outside the super- 
critical core (r ~ 0.4 pc, marked by a 'star'), where the gas is 
magnetically supported, the profile flattens significantly and is 
almost uniform. For radii in the range ~ 10""^ — 10~^ pc the slope 
d\nBz/d\nr ~ —1.1 (see Fig. ^jp). As the temperature increases 
at smaller radii, the slope approaches ~ — 1.4 down to the location 
of the hydrostatic core boundary, inside which the slope quickly 
asymptotes to zero. This zero slope is due to effective Ohmic dis- 
sipation, which quickly erases any spatial variation in the magnetic 
field. The value of the magnetic field inside the hydrostatic core is 
~ 0.2 G, in excellent agreement with the protosolar magnetic field 
as derived from meteoritic data I Stacev. Lovering & Parrv. .196k 



iHemdon & Rowelll974l ; lLevvlll988|y 

A closer inspection of Figures |9^ and |9j) reveals that the 



component of the magnetic field actually has a local maximum just 
outside the hydrostatic core at r ~ 10""" pc (~ 2 AU), where 
the inwardly advected magnetic flux has piled up and formed a 
'magnetic wall.' This local concentration of magnetic flux out- 
side the hydrostatic core boundary results in a magnetic shock, the 
consequences of which are dis cussed in Sections 13.3.31 and 13.3.41 
iTassis & MouschoviasI ( l2005bl) followed the formation and evolu- 
tion of a series of these magnetic shocks and found that accretion 
onto the forming protostar occurs in a time-dependent, spasmodic 
fashion. 

Fig. |9}; shows the radial profile of the maximum strength of 
the r-component of the magnetic field, at whichever height z it oc- 
curs above the midplane, at twelve different times (ti, t2, . . ., t\2) 
as in Fig. [8] Br is generated as a result of field-line deformation 
during core contraction. Inside the inner flat region where the mag- 
netic field is almost uniform, i3,.,max declines rapidly to very small 
values, reaching zero at the origin. The rapid contraction outside 
the hydrostatic core induces significant field-line deformation, and 
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Br, max has a sharp local maximum at the boundary of the hydro- 
static core, where its value becomes comparable to, although still 
smaller than, that of Bz (Br/ Bz = 0.99 at the location of Br, max). 
It is the case, nevertheless, that Bz > Br at all locations and for all 
times. 

The ratio of the local Alfven speed and the local isothermal 
sound speed, va/C, is shown in Fig.|9}l. The behaviour of this ra- 
tio may be understood from the fact that va/C oc rin Out- 
side the magnetically supercritical core (whose radius is marked 
by a 'star') k ~ and so va/C oc Wn . This accounts for the 
large values of va/C in the cloud envelope. In the magnetically 
supercritical core during the isothermal phase of contraction (when 
K ~ 0.47 and 7cff — 1), the Alfven speed varies only slightly with 
the isothermal sound speed and is approximately a constant of order 
unity. The fact that va becomes comparable to (or smaller than) the 
sound speed inside the supercritical core provides an explanation 
for th e thermalisation of l inewidths observed in m olecular cloud 
cores teaudrv et d]|l98ll : iMvers & Bensonlll983h . In the theory 
of magnetic star formation, the material motions responsible for 
the observed li newidths are attributed to long-wavelength, standin g 
Alfven waves ( lMouschoviaslll987l : iMouschovias & Psaltis|[l99^ . 
with a rem arkable quantitative agree ment between theory and ob- 
servations jMouschovias et ^|2006[) . 

A local maximum in ua/C occurs just outside the magnetic 
wall at r ~ 10 AU. There is a substantial decrease in va/C in- 
side the hydrostatic core, where the gas is primarily thermally sup- 
ported, due to the increase in 7eff and decrease in k; i.e., the gas 
begins to evolve almost adiabatically and the magnetic field begins 
to decouple from the matter. 



3.3.3 Forces 

Fig. [To] exhibits the relative magnitude of thermal-pressure forces, 
magnetic-pressure and tension forces, and gravitational forces dur- 
ing the evolution of the forming protostar. Fig. [Tob shows the ra- 
tio of thermal-pressure and gravitational forces. The gravitational 
force dominates the thermal-pressure force everywhere (since the 
cloud is thermally supercritical), except in the hydrostatic core 
(r < 2 AU). Just outside the hydrostatic core, there is an abrupt in- 
crease in the thermal-pressure force which is ultimately responsible 
for the formation of a thermal shock there (see Section [3.3.41 ). At 
time tw, the ratio of the thermal-pressure and gravitational forces 
inside the hydrostatic core is greater than unity. This is because the 
collapsing core has overshot its equilibrium, which results in radial 
(spherical) pulsations (see Section [3.3.4| (. 

Fig. 1 10b exhibits the ratio of magnetic and gravitational 
forces. The magnetic force is an appreciable fraction of the gravi- 
tational force everywhere except inside the hydrostatic core (where 
effective Ohmic dissipation is at work — see Section [3.3.6| (. This is 
the case even inside the dynamically contracting supercritical core. 
Hence, the dynamical contraction of the magnetically supercriti- 
cal core is significantly slower than free fall. The magnetic force 
has a local maximum with respect to the gravitational force at 
r « 10 AU, just outside the location of the maximum of the mag- 
netic field seen in Fig. [9^. There the abrupt increase in the magnetic 
force by a factor > 2.5 represents a magnetic wall and results in the 
formation of a magnetic shock (see Section [3.3.4b . 

In order to determine which agent is mainly responsible for 
diluting the effects of gravity at any given radius, we also give 
in Fig. llOb the ratio of thermal-pressure and magnetic forces. 
The magnetic force is greater than the thermal-pressure force (and 
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Figure 10. Radial midplane profiles of (a) the ratio of thermal-pressure and 
gravitational forces, (b) the ratio of magnetic and gravitational forces, and 
(c) the ratio of thermal-pressure and magnetic forces at different times, as in 
Fig. [8] The inner (outer) 'star' on a curve, present only after a supercritical 
core forms, marks the initial (final) radius of the supercritical core. 
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Figure 11. Radial midplane (solid line) and vertical symmetry-axis (dashed 
line) profiles of the velocity of the neutrals (a) in km s^^ and (b) in units 
of the instantaneous local sound speed, C. Different curves correspond to 
different times, as in Fig. [H The inner (outer) 'star' on a curve, present 
only after a supercritical core forms, marks the initial (final) radius of the 
supercritical core. 

hence provides the main opposition to gravity) everywhere except 
in the hydrostatic core, where the thermal-pressure force dominates 
the magnetic force by a few orders of magnitude. There is a local 
minimum in the thermal-to-magnetic force ratio outside the hydro- 
static core, where the magnetic field strength increases briefly. 

33.4 Velocities 

The velocity of the neutrals is given in Fig. [TT](a) in km s~" and 
(b) in units of the instantaneous local isothermal sound speed, C. 
For reference, the isothermal sound speed at a temperature of 10 K 
is 0.188 kms^". The ambipolar-diffusion-controUed, quasistatic 
phase is marked by small radial velocities (~ 0.01 — 0.03 km s""). 
Once a magnetically and thermally supercritical core forms, dy- 
namical contraction ensues inside the core and the radial velocity 
increases to ~ 0.5C ~ l.OC. It is noteworthy that the radial veloc- 



ity is ~ 0.1 km s " at r ~ 0.1 pc, in excellent agr e emen t with 
the observed infall motions of LI 544 bv lTafalla et"S] ( Il998l) . The 
predicted radial velocities also compare well with the observed in- 
fall motions near the starless molecular cloud core MC27 in Tau- 
rus, where speeds of 0.2 — 0.3 km have been detected at 
r ~ 2000 - 3000 AU (~ 10"^ pc) jOnishi, Mizuno & Fukuil 
1 19991) . The radial velocities outside the supercritical core, where 
the gas is magnetically supported, remain subsonic by a factor of a 
few. 

At approximately the same time that isothermality breaks 
down in the supercritical core, the radial velocities become super- 
sonic. The mass accumulating in the hydrostatic core causes an 
accelerated infall and the radial velocity approaches and exceeds 
~ 1 kms~". A shock forms near the boundary of the hydro- 
static core due to rapid increases in the local magnetic and thermal- 
pressure forces, and the velocity rapidly falls to zero toward the 
centre of the hydrostatic core. The large vertical velocities through- 
out the evolution result from the cloud collapsing along field lines 
and constantly responding to changes in the thermal-pressure gra- 
dient along magnetic field lines. At the boundary of the hydrostatic 
core, shocks appear along the (vertical) symmetry axis. Vertical ve- 
locities become positive after the core overshoots hydrostatic equi- 
librium, rebounds, and oscillates. Radial velocities remain negative, 
but form a thermal and magn etic shock outside the hydrostat ic core. 

The parameter study bv lTassis & Mouschoviaj j2007ch found 
that, if the onset of adiabaticity is delayed beyond its oft-assumed 
critical density lO"" cm^'', as is naturally the case in the RMHD 
simulation presented here, the maximum infall velocity of both the 
neutrals and magnetic field lines increases somewhat as the den- 
sity at which adiabaticity sets in increases. This is because the de- 
lay of the onset of adiabaticity implies accelerated infall for longer 
times. Using the results fro m a simulation in which isother mality 
was assumed to always hold, lTassis & Mouschoviaj ( 12007 d) found 
that the longer the gas remains isothermal, the weaker the mag- 
netic shock is. In other words, the enhanced 'pile-up' of matter 
and magnetic flux outside the hydrostatic core becomes increas- 
ingly absent. Given that the results of our RMHD simulation must 
necessarily be bracketed by those in which isothermal or adiabatic 
equations of state are assumed, the fact that our magnetic wall is 
more pronounced than found in their isothermal control run, yet 
less pronounced than found in their adiabatic run, is not surprising. 

Fig- shows radial profiles of the velocities and attachment 
parameters for the grains and the ions. (The electrons remain well- 
coupled to the magnetic field lines throughout the evolution.) The 
attachment parameter for species s is defined by 



and quantifies the degree to which species s is coupled to the 
magnetic field both directly (by gyrations of individual particles 
about field lines) and indirectly (through electrostatic attraction by 
oppositely-charged particles that are at least partially directly at- 
tached to the field lines). If As « 1, then species s is attached to 
the field lines, whereas if As -C 1, then species s is detached from 
the field lines and its motion follows that of the neutrals. 

Outside the magnetically supercritical core, all charged 
species and the neutral grains remain well attached to the magnetic 
field lines. It is only inside the supercritical core that one after an- 
other the charged species begin to detach from the field lines. The 
negatively-charged grains are detached from the field lines by the 
time te. the neutral grains by It, and the positively-charged grains 
by fg. The ions do not begin to detach from the field lines until t^. 
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Figure 12. Radial midplane profiles of (top row) velocities and (bottom row) attachment parameters for (left to right) positively-charged grains, neutral grains, 
negatively-charged grains, and ions. Different curves correspond to different times, as in Fig. 1] The inner (outer) 'star' on a curve, present only after a 
supercritical core forms, marks the initial (final) radius of the supercritical core. 



well into the dynamical contraction phase, and do not fully detach 
until iio. Indeed, at late times, the radial velocity profiles of these 
species shown in Fig.[T2]appear identical to those of the neutrals in 
Fig. [nil- 

A closer inspection of the attachment parameter profiles re- 
veals more detailed and interesting behaviour than one would ex- 
pect based solely on magnetic and collisional forces. These param- 
eters do not simply asymptote directly to zero, but rather asymp- 
tote first to a small but nonzero value before finally settling at 
zero. For example, Ag asymptotes to « 0.05 before falling to 
zero. Thus, negatively-charged grains remain mildly attached to the 
field lines even when simpler estimates (e.g., based on the crite- 
rion a;grgn 1) would suggest otherwise. This is a result of elec- 
trostatic attraction between negatively-charged grains and electron- 
shielded ions, which have yet to fully decouple from the magnetic 
field. Neutral grains also couple to the magnetic field by inelas- 
tic charge-capture processes. Slightly more complex behaviour can 
be seen for the positively-charged grains, which are repelled from 
the field lines by electron-shielded ions. Even electrostatic attrac- 
tion between ions and (fully-attached) electrons keeps the ions very 
mildly attached for - iio, though this effect is much subtler than 
for the grains because of the extremely low electron abundance at 
these high densities. 



The current density (in the midplane) due to species s is given by 
js,0 = risqsVs,,/!- Absolute values have been taken; the cusps in the 
positive-grain curves are due to the current passing through zero as 
it changes sign. This occurs at the same time that the positive-grain 
conductivities change sign (see Figs.|7}; and|7jl) and at the same 
spatial location that the positive-grain attachment parameter passes 
through zero and becomes slightly negative (see Fig.ll2t. 

At early times the contributions of electrons and ions to the 
electric current density are comparable and are approximately three 
orders of magnitude greater than the contribution from negatively- 
charged grains. By tg, the increasing disparity between electron and 
ion abundances (due to electron depletion onto grains) causes the 
ion current to become greater than the electron current by a factor 
~ 100 (recall that, at this stage, Xi/xe ~ 125 — see Fig. [6J5). 
Finally, in the innermost region, the current density decays sharply 
due to the increasing effect of Ohmic dissipation (see Fig.llSt. As 
discussed in Section [T2l grains become the main current carriers 
at roughly the same time as Ohmic dissipation becomes important. 
Their contribution to the electric current density is more than three 
orders of magnitude greater than that of the ions and electrons. 



3.3.5 Contribution to Electric Current Density by Different 
Species 

Fig.[T3]exhibits the contribution of different charged species to the 
total electric current density (in the midplane) as a function of ra- 
dius at different times. The current density js,<p (in the midplane) 
carried by species s may be obtained from equation (CI) of Pa- 
per I and written in terms of the total current density (shown in 
Fig.ll4ll and various components of the resistivity tensor: 



3.3.6 Ambipolar Diffusion and Ohmic Dissipation 

The relative importance of ambipolar diffusion and Ohmic dissi- 
pation as agents of magnetic diffusion is demonstrated in Fig. [15] 
which shows radial midplane profiles of the ratio of the ambipolar- 
diffusion and Ohmic-dissipation resistivities. At low densities and 
large radii (and height), ambipolar diffusion dominates Ohmic dis- 
sipation. The two processes become equally important only at a 
density ~ 7 x 10^^ cm"'', while at higher densities Ohmic dissipa- 
tion operates on a timescale shorter than that of ambipolar diffusion 
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Figure 13. Radial profiles of the magnitude of the electric current density in the 0-direction carried by (a) the positively-charged grains, (b) the negatively- 
charged grains, (c) the ions, and (d) the electrons. Different curves correspond to different times, as in Fig. [8] The inner (outer) 'star' on a curve, present only 
after a supercritical core forms, marks the initial (final) radius of the supercritical core. 



(see also Fig.|7^.) The region in which this occurs coincides with 
the region in which the electric currents diminish (r < 2 AU). 

Note that, by the time Ohmic dissipation becomes important, 
all species except the electrons have detached from the field lines. 
Hence, magnetic decoupling is essentially complete. Calculations 
that arrived at the conclusion that Ohmic dissipation is the cause of 
magnetic decoupling (e.g., Nakano & Umebavashilll986al fbl) were 
based on a value of the e-H2 cros s section which is to o large by a 
factor of about 100 (see review bv lMouschovias|[l996l . § 2.1). 



3.3. 7 Mass, Magnetic Flux, and Mass Infall Rate 

In Fig. |16t i we show radial profiles of the cumulative mass within 
a radius r. Within the thermally and magnetically supercritical 
core, there are three different regions that exhibit characteristi- 
cally different (logarithmic) slopes of the mass spatial profile: (1) 
Inside the hydrostatic core, both the density and the temperature 
are essentially uniform (see Figs. [8^ and [8]:), because the size 



of this region is smaller than the instantaneous value of the ther- 
mal critical lengthscale AT,cr(i) and therefore thermal-pressure 
forces rapidly smooth out any nonuniformity in these quantities. 
Consequently, dlnnn/dlnr ~ 0, from which it follows that 
d\n M/dlnr ~ 2 (see Fig. 1 16b, region r < 10^''' pc). (2) Just 
outside the boundary of the hydrostatic core, the infall velocity has 
its maximum value (see Fig. II lb , thereby tending to empty this 
region of mass and to cause a steepening of the density profile. 
Since relatively little mass exists in this region, the mass profile 
flattens (see Fig. 1 16b, region 10~® PC ^ 10~* pc). (3) Far- 
ther out in the supercritical core, where the geometry is disclike, 
we have that M = vrr^S oc p^"''^'^ , from which it follows that 
ainM/ainr ~ 2+(7cft /2)(ainnn/91nr) f» 0.9. By the end of 
the run, ~ 12 M© has accumulated within the magnetically super- 
critical core and ~ 0.006 M0 within the hydrostatic core. Hence, 
the often made statement that ambipolar diffusion can only form 
low-mass cores ( and low-mass stars ) is not correct. 

Similar behaviour is evident in Fig. 1 16b . which shows ra- 
dial profiles of the cumulative magnetic flux within a radius r. 
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Note, however, that the clear break in the mass-radius profile at the 
boundary of the hydrostatic core is hardly evident in the magnetic- 
flux-radius profile. The magnetic flux does not make the full tran- 
sition to spherical geometry like the mass does, since the magnetic 
field is decoupled from every species except the tenuous electron 
fluid (see Section D.3.4| l. The magnetic flux threading the hydro- 
static core is « 5 X 10"^ nG pc^ 4.8 x 10^** Wb. It is in- 
triguing that this is close to the magnetic flux of a typical Ap star 
such as 6 Aurigae, which has a radius 4.5 R© and magnetic field 
strength _B ~ 1 kG, giving a magnetic flux ~ 3 x 10^* Wb 
jvan Rensbergen. Hensberge & Adelmanlll984h . Ap stars are be- 
lieved to avoid the convective stage as protostars and so retain the 
magnetic flux leftover from their formation. Coupled with 9 Au- 
rigae's young age {r^ 200 Myr), this suggests that its ^ 1 kG 
magnetic field is a fossil field. 

The radial profile of the mass-to-flux ratio (normalised to the 
central critical value for collapse) is given in Fig.|16t;. A visual ex- 
amination of the profile at the last time output suggests that the 
curve exhibits a three-slope power law with a flat inner region, 
which corresponds to the flat-density uniform-magnetic-field core. 
However, a plot of the profile of the slope itself shows that there 
is no region in which it is constant (see Fig. I16tl). Conclusions 
based on visual inspections, especially of log-log plots, can be very 
misleading. Relatively abrupt changes in the slope separate regions 
governed by different physics. The first steepening (at small r) oc- 
curs just outside the hydrostatic core and marks the transition from 
the magnetically-decoupled hydrostatic core to the magnetically 
supercritical core, in which the magnetic flux is almost frozen in 
the matter. From radii ~ 10"* pc to a few xlO~'^ pc, the slope 
varies relatively slowly. In this region, the ions are still attached to 
the field lines and the mass-to-flux ratio hardly changes in time. 
Beyond the supercritical core, where magnetic forces continue to 
support the cloud envelope, the slope steepens again. 

The neutral density (Fig. [8ji) and velocity (Fig. II It profiles 
are used to calculate the mass infall rate, dM/dt, through any 
(cylindrical) radius r. The infall rate is given in Fig. 117b in units 
of Mq yr~\ It varies considerably both spatially and tempo- 
rally, a behaviour also noted in magnetic calculations by, e.g ., 
Mouschovias & Morton ( 1992a, b), Ciolek & MouschoviasI (1994), 
Basu (1997), and De sch & Mouschovias (200^7 and in nonmag- 
netic calculations by [Foster & Chevalieil ( Il993h . This behaviour 
contrasts sharply with that of the singular isothermal sphere model, 
which would pre dict a cons tant and uniform mass infall rate equal 
to 0.975 C^/G ( Ishulll977h . (The inclusion of magnetic fields to 
this model, which results in isothermal toroids, gives a greater 
but still time-in dependent and spatially uniform infall rate — see 
lLi&Shij|l996l .) 

In order to compare the mass infall rate more quantitatively 
with the predictions by the various self-similar solutions men- 
tioned earlier, we also plot in Fig. 117b the mass infall rate nor- 
malised to C'^/G, where C is the instantaneous local isother- 
mal sound speed. The horizontal dotted lines denote the con- 
stant mass infall rates predicted by th e Shu (1977) and Larson- 
Penston dLarsonl Il969l : IPenstonI Il969h solutions. At the bound- 
ary of the magnetically supercritical core the mass infall rate is 
initially 1.54 x 10"'^ M© yr"^ (= 0.96 /G), increasing to 
2.26 X 10"*^ Mq yr-i (= 1.41 G^ /G) by the end of the run. The 
maximum infall rate is 3.15 x 10~* Mq yr^^ (at r — 1.60 AU), 
whereas the maximum value of [dM / dt) / (C'* /G) = 27.08 (at 
r = 52.2 AU). That the early ambipolar-diffusion-controlled 
phase is characterised by a maximum mass in fall rate co mparable 
to the constant mass infall rate predicted by the lShtj il977h solution 
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Figure 14. Radial profiles of the magnitude of the total electric current 
density at different times, as in Fig. [8] The inner (outer) 'star' on a curve, 
present only after a supercritical core forms, marks the initial (final) radius 
of the supercritical core. 
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Figure 15. Radial profiles of the ratio of Ohmic and ambipolar resistivi- 
ties. The dotted line denotes the boundary when ambipolar diffusion and 
Ohmic dissipation become equally important magnetic diffusion mecha- 
nisms. Above (below) this line, Ohmic dissipation (ambipolai' diffusion) 
dominates. Different curves correspond to different times, as in Fig.[8] 

intuitively makes sense; the [shi] ( fT977h solution was obtained by 
assuming quasi-static (i.e., negligible velocity) initial conditions, 
which is similar to the quasistatic (i.e., negligible acceleration) evo- 
lution of the early, ambipolar-diffusion-controlled phase that exists 
before the creation of a magnetically supercritical core. Likewise, 
the mass infall rate in the Larson-Penston (1969) solution was ob- 
tained by assuming highly dynamical initial conditions, and so it 
is not surprising that the maximum mass infall rate found in our 
simulation, which occurs once a 'point mass' (the hydrostatic core) 
is formed in the central region, is bounded from above by the (spa- 
tially uniform) mass infall rate found in the Larson-Penston (1969) 
solution. 
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Figure 16. Radial profiles of (a) the cumulative mass within a radius r, (b) the cumulative magnetic flux within a radius r, and (c) the mass-to-flux ratio 
(normalised to the central critical value for collapse) at different times, as in Fig. [8] (d) Radial derivative of the mass-to-flux ratio, d\Ti{M /^-Q)/8\nr, at 
ti2 - The inner (outer) 'star' on a curve, present only after a supercritical core forms, marks the initial (final) radius of the supercritical core. 



3.4 Radiation Flux and Luminosity 

Under the FLD approximation, the radiation flux is given by 



(6) 



where Afld is the flux limiter, xr is the Rosseland mean extinction 
coefficient, and £ is the (frequency-integrated) radiation energy 
density (see Paper I). In Fig.[T8] we give the r-component of the 
radiation flux along the midplane (solid line) and the 2-component 
of the radiation flux along the symmetry axis (dashed line). To ease 
conversion to luminosity, 



L(^) 



(V) 



in the figure we have multiplied the flux by in and given it in 
units of Lq AU~^. (To obtain the luminosity at the hydrostatic 
core boundary in Lq, one simply multiplies the value of the plotted 
quantity evaluated at r ~ 2 AU by ~ 2^.) The radiation flux is sub- 
stantially greater along the symmetry axis than along the midplane, 



because there is much less absorbing material in the former direc- 
tion than in the latter. Near the boundary of the hydrostatic core the 
radial flux attains its maximum value ~ 2.5 x 10^" Lq AU~^ 
corresponding to a luminosity ~ 10^'^ Lq. 



4 SUMMARY AND DISCUSSION 

As shown by earlier calculations of our group, ambipolar diffusion 
in an initially magnetically subcritical cloud allows the neutrals 
to contract quasistatically due to their self-gravity through almost 
stationary magnetic field lines, redistributing mass in the central 
flux tubes of the cloud. Once the central mass-to-flux ratio exceeds 
the critical value for collapse (at a central neutral number density 
~ 10^ cm^'' under the conditions investigated here), a core begins 
to contract dynamically and is referred to as a supercritical core. 
Its contraction is dynamic, but slower than free fall, with its mag- 
netic flux almost frozen in the matter. Resurrection of ambipolar 
diffusion (at a density ^ 10^^ cm^'^) marks the onset of magnetic 
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Figure 17. Radial profiles of the mass infall rate through any (cylindrical) 
radius r in units of (a) Mq yr~^ and (b) the instantaneous local value 
of /G. Different curves correspond to different times, as in Fig. [8] The 
inner (outer) 'star' on a radial profile curve, present only after a supercritical 
core forms, marks the initial (final) radius of the supercritical core. The 
ho rizontal dot ted lines denote the constant mass infall rates predicted by 
the lShiJ fToTTh and Larson-Penston (1969) solutions. 



decoupling. All species except the electrons gradually, over several 
orders of magnitude density enhancement, detach from the field 
lines. Dust grains play a crucial role in the evolution of both the 
cloud and the core, and by the time these high densities are reached, 
the grains carry most of the electric charge, despite the fact that 
most grains are electrically neutral. At densities > 3 x lO"^ cm"'', 
the grains are also the dominant electric current carriers. These den- 
sities are comparable to those at which Ohmic dissipation domi- 
nates ambipolar diffusion (> 7 x lO"^ cm"''), which is found to 
occur in a central region of size ~ 2 AU. 

Once magnetic decoupling occurs in the central region, the 
magnetic flux in this region does not increase in time. Instead, 
infalling matter entering this region leaves its magnetic flux be- 
hind, thus creating a magnetic wall (at a radius i?waii ~ 10 AU), 
the boundary of the decoupling region. Inside, the magnetic field 



Figure 18. Radial midplane (solid line) and vertical symmetry-axis (dashed 
line) profiles of the radiation flux at different times, as in Fig. [8] To ease 
conversion to luminosity, we have multiplied the flux by in and have given 
it in units of Ln AU~^. 



becomes and remains spatially uniform, with magnitude Bdcc ~ 
0.2 G. It therefore follows that the protostellar/protoplanetary disc 
which will eventually form will be exposed to this Bdcc ~ 0.2 G. 
Hence, it may be no accident that measurements of the remanent 
magnetisation in meteorites imply magnetic fields in the early solar 
nebula of strength 0.1 - 0.2 G. 

Magnetic fields are also responsible for quickly establishing 
and maintaining a disclike geometry all the way down to radial 
lengthscales on the order of 10 AU, inside which thermal-pressure 
effects become important. This provides a natural explanation of 
observations of star-forming cores, such as L155 1-IRS5, that reveal 
a disclike geometry over a broad range of lengthscales. 

Throughout the entire collapse phase, the released gravita- 
tional potential energy escapes only as rapidly or slowly as allowed 
by the dust grains, which are mainly responsible for the opac- 
ity at temperatures < 1500 K. Although the gas remains strictly 
isothermal only for central densities < 10^ cin"^, the tempera- 
ture does not exhibit an appreciably rapid increase until densities 
~ 10^^ cm"'^ are attained in the core. A central temperature of 
100 K is reached when the central density is ~ 6 x 10^^ cm"^. 
For the densities considered here (< lO"* cm~^), the gas never 
evolves adiabatically. 

The dramatic temperature increase inside ~ 4 AU leads to the 
formation of a hydrostatic core, in which approximate balance is 
achieved between thermal-pressure and gravitational forces. At the 
boundary of the hydrostatic core (~ 2 AU), there are prominent 
shocks, particularly along the vertical symmetry axis where veloci- 
ties rapidly decrease from ~ — 1.2 km s~" to zero in only ~ 8 AU. 
The infalling matter overshoots the hydrostatic equilibrium state 
and subsequently oscillates about it. A simple extrapolation of tem- 
perature based on our results indicates that a temperature of 1000 K 
will not be reached until a central density ~ 3 x 10^^ cm"'^. There- 
fore, the recoupling of the gas to the magnetic field, via thermal 
ionisation and grain sublimation, will most likely not occur until 
central densities of at least several x lO"^ cm"^ are attained. 

At the end of the calculation, the mass and magnetic flux in 
the hydrostatic core are 0.006 M0 and 5 x 10~^ ^G pc^, respec- 
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lively. The mass-to-flux ratio in the central flux tube is ~ 80 times 
the critical central value for collapse (or ~ 169 times the uniform 
thin-disc critical mass-to-flux ratio). The luminosity at the hydro- 
static core boundary is ~ 10""^ L0. The mass infall rate is highly 
nonhomologous and time-dependent, rising from ~ C'-^/G ~ 
1.5 X lO"'^ M© yr~^ at the initial boundary of the magnetically 
supercritical core to a maximum value 3 x 10"'' M© yr~' 
at a radius of 1.6 AU, near the boundary of the hydrostatic core. 
The maximum value of the mass infall rate normalised to C'^/G is 
~ 27. The quantitative similarity between this value and the (spa- 
tially constant) mass infall rate associated with the Larson-Penston 
self-similar hydrodynamic collapse solution is no coincidence. By 
the time a spherical hydrostatic core is formed in our simulations, 
the infall velocities are dynamic and the magnetic field in the core 
has effectively decoupled from the matter. 

Despite a concerted effort to treat the nonisothermal stage of 
magnetic star formation as rigorously as possible, there still re- 
mains a great deal of work to be done, including both improvements 
to and extensions of what has already been accomplished here. 

First, dust coagulation ought to be computed in a time- 
dependent fashion. Unfortunately, the exact details of how and 
even when grain growth occurs in protostellar cores and discs are 
poorly constrained both experimentally and observationally. Due 
to recent improvements in instrumentation, however, observational 
prospects are becoming increasingly auspicious. A recent compre- 
hensive study of the circumstellar disc surrounding the T Tauri star 
IM Lupi using photometry, spectroscopy, millimetre interferome- 
try, and multi-wavelength imaging has obtained qua ntitative evi- 
dence of dust evolution in the disc i lPinte et alj|2003) . Similar re- 
sults have also been found for a sample of protoplanetary discs 
around 21 Class II T Tauri s tars in the Taurus- Auriga star form- 
ing region jRicci etal.ll2009h . The details, nevertheless, may not 
be very important, since the effect of grain size on the tempera- 
ture evolution is expected to be minimal, at least for the range of 
densities studied here; different grain sizes share roughly the same 
opacity for temperatures less than ~ 400 K. On the other hand, the 
strength of th e magnetic field at decoupling is a sensitive function 
of grain size jDesch & Mouschoviasll20()lr) . One possible way to 
circumvent the uncertainties of grain growth, at least at the present 
time, would be to adopt a nonuniform grain size distribution (such 
as the MRN distribution) and rerun the numerical simulation pre- 
sented here. As explained in Paper I, we have already modified the 
Zeus-MP code to handle such a distribution. Results will be pre- 
sented in a separate paper, as part of a parameter study. 

Second, although rotation was justifiably ignored for the den- 
sities investigated here (see Paper I), it cannot be ignored in a study 
of the more advanced stages of evolution. It is essential for under- 
standing the formation and evolution of protoplanetary discs. It is 
unknown at present whether magnetic braking plays a role in the 
formation of such discs and in their angular velocity structure or 
whether turbulent transfer of angular momentum may be impor- 
tant. 

A natural extension as well as an important test of our present 
work and predictions is the production of synthetic observations to 
be compared with ongoing actual 3 mm dust and N2H+ observa- 
tions of the very young Class source HH211-mm by colleagues 
using the Combined Array for Research in Millimeter Astronomy 
(CARMA). Such close interplay between theory and observations 
is now under way. 
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